Ultracold bosons with the synthetic three-dimensional spin-orbit coupling in an 

optical lattice 
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We study ultracold bosonic atoms with the synthetic three-dimensional spin-orbit (SO) coupling 
in a cubic optical lattice. In the superfluidity phase, the lowest energy band exhibits one, two or four 
pairs of degenerate single-particle ground states depending on the SO-coupling strengths, which can 
give rise to the condensate states with spin-stripes for the weak atomic interactions. In the deep 
Mott-insulator regime, the effective spin Hamiltonian of the system combines three-dimensional 
Heisenberg exchange interactions, anisotropy interactions and Dzyaloshinskii-Moriya interactions. 
Based on Monte Carlo simulations, we numerically demonstrate that the resulting Hamiltonian 
with an additional Zeeman field has a rich phase diagram with spiral, stripe, vortex crystal, and 
especially Skyrmion crystal spin-textures in each xy-p\ane layer. The obtained Skyrmion crystals 
can be tunable with square and hexagonal symmetries in a columnar manner along the z axis, and 
moreover are stable against the inter-layer spin-spin interactions in a large parameter region. 

PACS numbers: 67.85.-d, 05.30.Jp, 71.70.Ej, 37.10.Jk 



I. INTRODUCTION 

Spin-orbit (SO) coupling plays an important role in 
condensed matter physics, especially in the newly discov- 
ered quantum spin Hall effects and the related topological 
orders [J ■ Recent theoretical proposals 0, H| and ex- 
perimental realization [BHU of non-Abelian gauge fields in 
ultracold atoms with the optical dressing technique open 
another door to explore SO coupling in controllable sys- 
tems. The bulk gases of weakly interacting bosons with 
the synthetic two-dimensional (2D) Rashba SO coupling 
in homogenous cases and in trapping potentials have been 
widely studied and theoretically shown to exhibit exotic 
many-body ground states [9j, some of which have no di- 
rect analog in solid-state systems [l0rll3| . For example, 
an SO-coupled spinor condensate will spontaneously de- 
velop a plane wave phase or spin-stripe structure depend- 
ing on the weak interaction energy [l(| [H| ; and in the 
presence of strong trapping potentials, it will exhibit half- 
quantum vortex states and Skyrmion patterns pd| - [l3j . 

Recently, physics of a SO-coupled bosonic gas loaded in 
a 2 D op tical lattice (OL) has attracted considerable inter- 
est |l4j--23] . This system can be described by an extended 
two-component Bose-Hubbard (BH) model [24], [25[, in 
which the SO coupling can significantly affect the quan- 
tum phase transition from a superfluid to a Mott insu- 
lator (MI) [HElHii. More interestingly, the effective 
spin Hamiltonian of the system in the deep MI regime 
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contains the so-called Dzyaloshinskii-Moriya (DM) inter- 
action term [26[ , which comes from the SO coupling and 
may lead to some novel magnetic phases |16l - fl9l |. such 
as Skyrmion crystals [l6j . However, the corresponding 
three-dimensional (3D) system is yet to be explored. 

On the other hand, the 3D analog of SO coupling in 
cold atoms has been proposed to be experimentally re- 
alized by using optical dressing schemes [13, [H| and by 
exploiting laser-assisted tunneling (29[ in OLs [3(J [Hf . 
The 3D SO-couplings are less explored in contrast to the 
standard 2D Rashba and Drcsshauls ones in the context 
of condensed matter physics [271, but is now attract- 
ing more and more interests |2|. l32l |33| for investigat- 
ing 3D topological insulators Ll topological superfluid- 
ity (33[, and Weyl semi-metals (32J. Very recently, sev- 
eral pieces of theoretical work study the ground state of 
a weakly-interacting two-component Bose-Einstein con- 
densate (BEC) with the 3D SO-coupling in the contin- 
uum [34 36] . It is shown that the density distribution 
of the ground state BEC can also exhibit the interest- 
ing Skyrmion structure, which is moreover a 3D counter- 
part characterized by a 3D topological winding number 
[341 135] ]. Some other schemes have also been proposed 
to create Skyrmions in multicomponcnt BECs in the ab- 
sence of synthetic SO-couplings and OLs [37j ■ So it would 
be worthwhile to search the stable Skyrmion (crystals) in 
an OL system with 3D SO-coupled bosons. 

In this paper, we investigate ultracold bosons with 
the synthetic 3D SO coupling in a square OL. We first 
look into the weakly interacting superfluidity case. In 
this case, the lowest energy band exhibits one, two or 
four pairs of degenerate single-particle ground states re- 
lated to the SO coupling, and each pair contains op- 
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posite wave vectors with values depending on the SO- 
coupling strengths. This can give rise to the conden- 
sate states with spin-stripes for the weak atomic inter- 
actions. We then focus on the deep MI regime with one 
atom per lattice and derive the effective spin-spin inter- 
action Hamiltonian of the system. The spin Hamilto- 
nian is a combination of three-dimensional Heisenberg 
exchange interactions, anisotropy interactions and DM 
interactions. Based on Monte Carlo simulations, we nu- 
merically demonstrate that the resulting Hamiltonian 
with an additional Zccman field has a rich classical phase 
diagram with spiral, stripe, vortex crystal, and espe- 
cially Skyrmion crystal spin-textures in each xy-plane 
layer. We find that the obtained Skyrmion crystals can 
be square or hexagonal symmetries with experimentally 
tunable parameters by varying laser-atom interactions. 
Moreover, the Skyrmion crystals in a columnar manner 
of extending along z axis are stable against the inter-layer 
spin interactions within a large parameter region. This 
cold atom system with high controllability in the effec- 
tive spin interactions may provide an ideal platform to 
further study exotic quantum spin models and find new 
phases of matter. 

The paper is organized as follows. The next section 
(Sec. II) introduces an extended BH model which de- 
scribes cold bosons with the synthetic 3D SO coupling in 
a square OL. In Sec. Ill, we briefly analyze the single- 
particle energy band and the properties of the weakly 
interacting superfluidity phase. In Sec. IV, we derive the 
spin Hamiltonian of the system in the MI regime, and 
present its rich classical phase diagram with interesting 
spin configurations. In particular, we study the profiles 
and stability of the Skyrmion crystals in the system. A 
brief discussion and short conclusion are finally given in 
Sec. V. 



II. MODEL 

Let us consider an atomic gas of pseudospin-1/2 ultra- 
cold bosons loaded into a 3D cubic optical lattice with the 
synthetic SO coupling. The single-particle Hamiltonian 
of the system is written as 



Ho = — + n x a x p x + KydyPy + K z a z p z + V(x, y, z), (1) 

where m is the atomic mass, p is the momentum op- 
erator, Krj with rj = x,y,z is respectively the strength 
of SO coupling along the rj axis, and a x y z are the three 
Pauli matrices. Here the cubic optical lattice V(x, y, z) = 
V v sin 2 (k ri) is formed by three standing- wave laser 
beams with the same wave number k . Thus the lattice 
spacing is a = n/ko. Hamiltonian (fT]) can be rewritten 
as Hq = (p — A) 2 /2m + V up to a constant, where the 



non-Abelian gauge potential A = —m(K x <r x , k ,,o 
This is corresponding to a 3D SO coupling 

We consider the system in the tight-bindin g re gime 
which is reachable in realistic experiments 24, 25]. Un- 



der this condition, in the presence of such a non-Abelian 
gauge field, the bosons can be described by an extended 
single-band BH Hamiltonian in terms of Peierls substitu- 
tion EE [13: 
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where al (a- ha ) creates (annihilates) a spin-er (a =t,i) 
boson at the site i. The first term in Hamiltonian ([2]) 
describes atomic hopping between neighbor lattice sites, 
with tf) representing the overall hopping amplitude in the 
absence of the synthetic SO coupling. The 2x2 matrix 
lZ ri = exjp(—j-A n a) is the Peierls substitution along di- 
rection fj with respect to the gauge potential. We rewrite 
IZn = exp(i0 I) cr ?? ) = cos#, ; l + ia n sin^ with dimension- 
less SO-coupling strength 9 n — TrmK v /hko. The diagonal 
and off-diagonal terms in the matrix respectively refer to 
the spin-conserving hopping and spin-flip hopping due 
to the SO coupling in the xy plane. The second term 
in Hamiltonian ^ denotes the atomic repulsive interac- 
tions, which is given by 

Vint = \ E U TT'a\^a\ t(7 ,ai^,ai t(J , (3) 



where U aa > is interaction strength between spins a and 
a' ' . The atomic interactions are almost spin-independent 
in experiments in the absence of Feshbach resonances 
[E IH, and hence we assume = = U and 

Ufi = E7j,-f = all with awl. In fact, the slight dif- 
ference between the intraspecies and interspecies inter- 
action strengths will help to select the degenerate many- 
body ground states of the system in the weak interacting 
superfluidity phase. 

This extend 3D BH model also exhibits the superfluid- 
ity and the MI phases for the weak and strong atomic in- 
teractions compared with the hopping energy 24, 25] , re- 
spectively. The quantum phase transitions between them 
are affected by the SO coupling in a similar manner as 
that in the 2D cases [IE [IE HO]- Thus we just consider 
the system in the two interaction limits in this work, and 
focus on the effects of synthetic 3D SO-coupling in the 
superfluidity and the MI phases. 



III. SUPERFLUIDITY STATES 

In this section, we consider the weakly-interacting su- 
perfluidity phase. In this regime, the hopping term dom- 
inates in Hamiltonian @. We first look into the hopping 
Hamiltonian Ht =% — V- ln t to obtain the energy band of 
the system, and then briefly discuss the effects of weak 
atomic interactions. The corresponding Hamiltonian in 
the momentum space can be written as 
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By using spatial Fourier transformations on Ht, we 
can obtain = H^ x + H^y + Hf. z with H^„ — 

— 2trj cos 9 V cos(k n a)l + 2t v s'm9 n sm(k v a)a v . Diagonal- 
izing ifk yields the energy structure of the system in the 
vanishing interaction limit: 

E^ = — 2 t v cos 9 V cos(k v a) 

±2 i^sm 2 9 v sm 2 (k ri a). (5) 

The lowest energy states in the resulting lower Bloch 
band E^ present candidates for the many-body ground- 
state of the bosonic gas with weak interatomic interac- 
tions. The Bloch momentum of these states denoted 
by k = (fcQ,fco7^o) can be directly obtained by solv- 
ing the equation (E^) |k=k = for minimizing E^ 
with the specific parameters t v and 9^. We find that 
there are possibly one, two or four pairs of degenerate 
minima in the energy band for different 9 V in the non- 
vanishing SO coupling cases. Each pair contains op- 
posite wave vectors with values depending on 9 n . For 
example, when 9 X = 9 y = O.lir and 9 Z = OAn, two 
degenerate minima locate at ko = (0, 0, ±9 z /a); when 
9 X = 9 y = 0.47T and 9 Z = O.lir, four degenerate min- 
ima locate at (±£i/a, ±£i/a, 0) with tan£i = tan 9 X /V2; 
when 9 X = 9 y = 9 Z = OAn, eight degenerate minima lo- 
cate at (±£2/0, ±£2/0, ±£2/0) with tan £2 = tan 9 X /\^3. 
This is in sharp contrast to the continuum case where the 
rotationally symmetric dispersion has an infinite ground 
sates degeneracy forming an SO sphere [13] ■ Owing to 
the reduction of degeneracy, the Bose condensates with 
the synthetic 3D SO-coupling in an OL would be more ro- 
bust against quantum fluctuations than their bulk coun- 
terparts. 

For the weakly interacting cases (i.e., U <C t n ) within 
the Gross-Pitaevskill (GP) approximation, the analy- 
sis of the ground state (condensate) wave-fucntion in 
this system is in parallel to those of the counterparts 
in 2D OL s as well as in 2D and 3D continuum 

cases [l(J, l42|. Hence we just present the conclusions 
here without detailed calculations. The condensate wave- 
function can be written as a superposition of all pos- 
sible single-particle (plane-wave) wave-functions of the 
the lowest Bloch states discussed above, and the cor- 
responding superposition coefficients are determined by 
minimizing the mean- field GP interaction (i.e., density- 
density interaction) energy [T(], EH- The GP interac- 
tion energy can be divided into the spin-independent and 
spin-dependent parts. The spin-independent term yields 
the same Hartree-Fork energy for any different selection 
of superposition coefficients, but the spin-dependent one 
with respect to a selects the coefficients for minimizing 
itself. For a < 1 , only one of the lowest degenerate single- 
particle state (whose numbers can be two, four or eight) 
is occupied, and thus the condensate wave-function is a 
plane-wave state with a finite momentum. On the other 



hand, for a > 1, one of the paired degenerate states 
(whose numbers can be one, two or four) are occupied 
with equal superposition coefficients, giving rise to the 
condensate states with spin-stripe density distribution 
[lol [Tl1 | . The structure of the spin-stripe is dependent 
on the vector ko and hence is tunable by the synthetic 
SO coupling. We note that these ground states are still 
degenerate except the case of stripe states with only two 
degenerate minima in the Bloch band. To further remove 
this accidental degeneracy, one should consider quantum 
fluctuations [43j . 
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FIG. 1: (Color online) Classical phase diagram of the reduced 
2D layer version of the spin Hamiltonian with t z — 0, ob- 
tained from Monte Carlo simulations. This is related to the 
bosons loaded in the OL in the deep MI regime. FM denotes a 
ferromagnet with magnetic moment along one direction (2 di- 
rection here), and correspondingly AFM an anti- ferromagnet 
with staggered moments. SP denotes the spiral phase which 
has a spiral wave in its spin configuration. ST denotes the 
stripe phase. SKX and VX respectively denote a 2D Skyrmion 
crystal and a vortex crystal. Same typical spin configurations 
of unusual phases (SP, ST, SKX, and VX) are shown in Fig. 
2 and Fig. 3. 



IV. SPIN MODEL IN MOTT-INSULATOR 
REGIME 

A. Effective spin Hamiltonian 

In this section, we turn to consider the system in the 
MI phase. We are interested in the MI regime with 
U 3> t-q and nearly unit atom per lattice site. In this 
case, the atoms are localized in individual lattices and 
the nearest-neighbor hopping can be treated as a per- 
turbation, leading to an effective spin Hamiltonian [44j j . 
To obtain the spin Hamiltonian of the system, one can 
begin with a two-site problem [45| . In the zero order, 
the system is described by the interaction Hamiltonian 
(J3j> , and the ground state manifold for the two-site prob- 
lem with one atom in each site composes four degenerate 
zero-energy states {| f;t),l t;!),l l;t),l Ul)}- Here 



4 



we have assumed uniform in-site energy and chosen it as 
the energy base. The exchange of two atoms in differ- 
ent sites does not require energy, and hence the single 
atom hopping should be eliminated in the second order 
with respect to the ratio t v /U. In this progress, there 
are six excited states {| U\fy, |0;tj), I tt;0), |0;tt>, I U 
; 0), |0; \.\.}} with an energy U. The hopping perturbation 
described by Ht couples the ground-state manifold and 
the excited-state one. The resulting effective Hamilto- 
nian up to the second order of perturbation reads Q |45[ 



(H cS ) 



71/ 



(G) 



where f3 and v label the four states in the ground-states 
manifold and 7 labels the six excited ones. 

After obtaining the two-site effective Hamiltonian from 
Eq. ([5]), it is straightforward to extend it to the lat- 
tice counterpart by introducing nearest neighbor hop- 
ping in the whole lattice. It is convenient to write the 
lattice effective Hamiltonian in terms of isospin opera- 
tors Si = (S?,S?,S?) with Sf = |(4t a U + »u ai 't)> 
S? = -|(aj jt ai,4.-a[j_ai,t), andSf = ^(Oi+Oi^-aJ^oy,). 
The resulting spin Hamiltonian of this 3D system in in 
the deep MI region is then given by 



H s 



i 

E 

i 

E 



Si ■ (jxSi+sb + JyS\+y + J Z S\- 

(^K x Sf S? +x + K y S\S\ + ~ + K Z S?S{ 



{d x §\ x S\-\-x ■ 
+ D z SiX Su 



DySi x Si +y ■ y 



where J, 

At 



At 2 

^cos(20„) 



K„ = -rf sin 2 



''t- 



is ^w\"VT}J> — u 

^-Bm{26 7l ) are spin interaction strengths. Hamiltonian 
(0) describes a generally anisotropic 3D spin-spin inter- 
action system. It combines the Heisenberg exchange in- 
teraction as the first term, the anisotropy interaction as 
the second term, and 3D DM spin interaction as the last 
term. Note that all the spin interaction strengths in these 
terms are dependent on the laser beams which generate 
the OL and the SO coupling, and hence they are tunable 
in experiments. 

To proceed further, we introduce an effective Zeeman 
term to Hamiltonian (O , leading to the total spin Hamil- 
tonian 



^E 5 i 



(8) 



This Zeeman term can be easily achieved by applying 
an additional external field to the pseudospin-1/2 atoms. 
For the pseudospin states that are usually two atomic hy- 
pcrfine states, the external field can be simply a real mag- 
netic field [45|, |4(| . If the pseudospin states are dressed 
states, one can use an appropriately designed laser field 



to generate it 47]. Thus h z is also a tunable parame- 
ter. We assume that the strength of Zeeman field here 
h z <^ U but is comparable with tlL/U. Under this con- 
dition, the parameters in Hamiltonian (J7J are approxi- 
mately unchanged. 



(a) 

16 
12 
>> 8 
4 






(M„ 



1.0 



I 



-0.5 



(e) 

M 



■ 



(d)„ 



4 8 12 16 
X 



11 iO 1 * ^ nJ^ j< ) 



X X X X X _X X X X 



X 



12 



16 



a) 


1 






e) 




d) 


• 
• 



-1 














1 
kjn 


- 


1 

k x /;i 


- 1 - 


1 

k x /7T 


- 1 - 


1 

k x /jt 



FIG. 2: (Color online) The spin configurations (distribution 
of S) in the xy-plane with z = 0: (a) SP {6 = 0.104, hz = 0); 
(b) VX (6» = l,h z = 0); (c) ST (0 = 1.24, h z = 0); (d) ST 
with complicated structure (6 — 1.56, h z — 2.34), which is not 
in the phase diagram, (e) The corresponding spin structure 
CO factors ( see the text) of the spin configurations in (a-d) 
and D = are f rom left to right, respectively. 



B. Numerical results from Monte Carlo simulations 

Below we explore the low-temperature phase diagrams 
and ground states of the effective spin Hamiltonian %J 
[see Eq. (|8J)], and we also intend to find stable Skyrmionic 
spin-textures via Monte Carlo simulations. In this clas- 
sical approximation we treat the spins Si as classical 
unit vectors and aim to find the spin configurations {Si} 
for minimizing the energy. For simplicity, we assume 
isotropic parameters in the xy plane, i.e., t x = t y = t 
and X = 8 y = 9, and take At 2 /U as energy unit hereafter, 
such that J x — J y = 008(26*), K x = K y = 2 sin 2 9 and 
D x = D y = sin(20). We also assume t z = X±t and 9 Z = 
X 2 9, such that J z = \\ cos(2A 2 6»), K z = 2\\ sin 2 (A 2 6») 
and D z = A 2 sin(2A26 ) ). Calculations were mostly carried 
out for an 18 x 18 x 18 lattice with periodic boundary con- 
ditions. Metropolis Monte Carlo algorithm [48[ was used 
throughout the calculations with 5 x 10 6 (and 5 x 10 5 for 
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an xy-layer with 18 x 18 lattice) sampling steps at each 
annealing process with fixed low temperature T = 0.005 
(in unit of At 2 /U). Some checks on each xy plane with 
sites 36 x 36 were performed to ensure consistency. The 
numerical results were also confirmed to be nearly the 
same for the open and periodic boundary conditions. 










X 


♦ 






■ 










as) * 




4 


I 



I 0.02 
-0.004 
-0.02 
-0.04 
-0.06 



FIG. 3: (Color online) Properties of the 2D Skyrmion crys- 
tals. Typical spin configurations in the Skyrmion crystal 
phase with (a) square and (b) hexagonal symmetry, and the 
color scale shows the magnitude of out-plane component S z . 
(c) and (d) are the corresponding spin structure factors of 
the spin configurations in (a) and (d), respectively, (e) and (f) 
are the local density x\ ( see the text) of the square and the 
hexagonal Skyrmion crystals, respectively. The parameters 
are 9 = 0.726 and h z = 0.910 for the square Skyrmion crys- 
tals; 6 — 0.592 and h z = 0.756 for the hexagonal Skyrmion 
crystals. 

We first consider the case of t z ~ in the Hamiltonian 
which can be realized by increasing the intensity of 
the lasers that generate the periodic lattices along z axis 
and freeze the atomic motions in this direction. In this 
case, the 3D system is equivalent to a collection of inde- 
pendent 2D plane layers along z axis, and thus the spin 
configurations in each layer are always the same in the 
ground states. We have confirmed this point in our nu- 
merical simulations. So we can first look into a single 
layer and figure out the phase diagram of Hamiltonian 
with t z = 0. In this limit, the s pin Hamiltonian is 
similar to those in the previous work [1(| [I?} without the 
additional Zeeman term. We have checked that in this 
case our results are consistent with those in Refs. [l6l.[l~7| 
when the parameter regions are overlapped, i.e., h z = 0, 
t x = t y , and 9 X = 9 y . 



We obtain the classical ground state phase diagram 
of the spin Hamiltonian with t z = as shown in 
Fig. 1. Apart from the conventional ferromagnetic and 
anti-ferromagnetic phases, there are four unconventional 
phases having interesting spin configurations in the xy- 
plane: spiral, stripe, vortex crystal and Skyrmion crys- 
tal. The presented phase diagram shows a rich interplay 
between different magnetic orders, and the parameter re- 
gion of the six phases can be found in Fig. 1. We should 
note that their boundaries between different phases may 
be unprecise in the level of quantum phases, however, 
this classical approach can be an efficient method to de- 
termine possible phases fTsj - TTsj . 

The spin textures in the four unconventional magnetic 
phases have non-trivial structures and can be character- 
ized by their spin structure factors S£ = | S^e lk " ri \ 2 
with = (Si,Sf,0). Figures 2(a)-(d) show some typ- 
ical spin configurations in the spiral, vortex crystal and 
stripe phases, and figure 2(e) shows their spin structure 
factors S^r in the momentum space (i.e., k x -k y plane), 
with the spots denoting the peaks of the spin structure 
factors. In Fig. 2(e), S£ exhibits a peak at (0,0.127r), 
corresponding to the spins spiralling along y axis with 
the wave number 0.127T as shown in Fig. 2(a). Simi- 
larly, the spins forming a vortex-crystal configuration in 
Fig. 2(b) has S£ peaks at (0.7tt, 0) and (0, 0.7tt), and the 
stripe spin configuration in Fig. 2(c) has a single peak 
at (it, 0), which means that the spins are staggered along 
x axis but are parallel along y axis. We also find that 
the ground states may exhibit more complicated stripe 
spin-textures outside the parameter region of the phase 
diagram, with an example being shown in Fig. 2(d). 

The Skyrmion crystal phase in the phase diagram, re- 
ferring to the phase where the spins presenting an array of 
2D Skyrmions (see Fig. 3), has a large parameter region 
when varying the SO-coupling strength and the Zeeman 
field h z . The array of Skyrmions is anisotropy in general, 
and can be square and hexagonal symmetric for proper 
SO-coupling strengths. Figures 3(a) and (b) show the 
typical spin configurations in the Skyrmion crystal phase 
with square and hexagonal structures, respectively. The 
corresponding spin structure factors 5*^ shown in Figs. 
3(c) and (d) also reflect their symmetries in the momen- 
tum space. The 2D Skyrmion crystals obtained in this 
model are distinguished from those in Ref. [l6| not only 
in their spin configurations but also in the mechanism. 
The 2D model there [l6[ contains no Zeeman field, how- 
ever, the Zeeman field is crucial here. As seen from the 
diagram, the Skyrmion crystals are generated with the 
help of a Zeeman field from the spiral waves (which are 
formed via the competition between the DM interaction 
and the Hciscnberg exchange interaction). Such kind of 
Skyrmion crystals have been explored by numerical cal- 
culations and in experiments in chiral magnet materials 
[li^ill, but the controllability in the materials [U EH 
is low in contrast to that in the cold atom system. This 
system is clean and the paraments, such as 9 and h Zl 
are widely tunable via adjusting laser-atom interations 
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[24 |25[ . So the density of Skyrmions and the symmetry 
of the Skyrmion crystal can be well controlled by varying 
the two parameters. For instance, increasing the Zeeman 
field can lead to the increase of the Skyrmion density 
at the beginning and the melting of the Skyrmions up 
to certain levels, and varying the SO-coupling strength 
can change the distribution of Skyrmions, from generally 
anisotropy to square or hexagonal symmetry as shown in 
Figs. 3(a) and (b). 

To further characterize the Skyrmion crystals, we in- 
troduce the local density of Skyrmions \\ at lattice site 
i in each xy-plane as (3S| 51] 



1 



Si ■ (S i+£ x + Si ■ (Si-i x , (9) 



which is the discretization counterpart of the well-known 
topological charge density S- (d x Sx d y S)/4:iT for the con- 
tinuum case [3J]. For a single localized 2D Skyrmion 
here, its topological winding number given by Wid — 
ccll Xi P ms the sign of its pole (i.e., here S z = —1 
at the Skyrmion cone as shown in Fig. 3) equals to unit 
in the continuum limit. Note that for an ordinary vor- 
tex, this topological number equals to zero. The winding 
number is stable with respect to the discretization, as 
for a lattice layer with size L x L, the fluctuation (er- 
ror) is on the order of 0(4n 2 /L 2 ) In Figs. 3(e) and 
(f), we show local density of the square and the hexago- 
nal Skyrmion crystals. The winding number of a single 
Skyrmion in the two cases is numerically computed to be 
nearly minus one. 

Finally in this section, we consider the inter-layer spin- 
spin interactions along z axis and check the stability of 
the square and hexagonal Skyrmion crystals obtained 
previously [see Fig. 3(a) and (b)] with respect to the 
parameters Ai and A2. Our numerical results are shown 
in Fig. 4. We find that the Skyrmion crystals with 
both square and hexagonal structures are stable and have 
a large region in the parameter space (Ai,A2). Since 
the parameters Ai and A2 are position-independent, the 
spin configurations still exhibit the same distribution of 
Skyrmions in each xy layer, which has been confirmed 
in our numerical calculations. That is to say, the 2D 
Skyrmion crystal states extend along the z axis in a 
columnar manner [501 ] . In the parameter space (see Fig. 
4), there is another phase, i.e. the ferromagnetic phase. 
This demonstrates that the layer Skyrmion crystals can 
be melt into the conventional ferromagnet by spin-spin 
interactions along z axis in Hamiltonian ©. We note 
that near the boundary between the two phases, the 
Skyrmions and ferromagnet coexist in the spin config- 
urations and hence the Skyrmion crystals are ill-defined 
there. However, in most of the region denoted by SKX in 
Figs. 4(a) and (b), the square and hexagonal Skyrmion 
crystals remain. In our numerical calculations done by 
the classical Monte Carlo method, we do not find a topo- 
logically non-trivial spin structure forming a genuine 3D 
Skyrmion or its corresponding crystal [511 ] for the spin 
Hamiltonian (J7J in a large parameter region. Whether 



such a ground state of 3D Skyrmion crystal can exhibit in 
the context of chiral magnetism is still an open question 
[39l [50l | . This cold atom system with widely tunable pa- 
rameters may provide a better platform for exploring the 
Skyrmion physics and searching for 3D Skyrmion crys- 
tal states, which would be an interesting challenge in our 
further studies. 




FIG. 4: Stability of layer Skyrmion crystals with respect to 
the spin-interaction parameters Ai and A2 along z axis. The 
parameter regions for (a) square and (b) hexagonal Skyrmion 
crystals denoted by SKX are both wide. Another region de- 
noted by FM is the ferromagnetic phase. The parameters are 
9 = 0.726 and h z = 0.910 in (a); 6 = 0.592 and h z = 0.756 in 
(b). 



V. DISCUSSION AND CONCLUSION 

Before concluding this paper, we briefly discuss some 
methods for detecting the superfluidity states and the 
spin configurations in the deep MI phase. The plane- wave 
and the strip superfluidity phases correspond to BECs at 
a single finite momentum and at a pair of opposite mo- 
menta, respectively. The values of the momenta depend 
on the SO-coupling strength as discussed in Sec. III. The 
standard time-of-flight imaging measurement can reveal 
the condensate peaks at the nontrivial momentum points 
[25I ] . which provides direct signatures of the two super- 
fluidity states. The different magnetic orders presented 
above can be detected by the optical Bragg scattering 
for atoms in OLs [Hj], as the peaks in the Bragg spec- 
troscopy directly reveal their spin structure factors. An- 
other way to measure the spin configurations is using the 
spin-resolved in-situ imaging technique [Hj]. The idea is 
to implement a high-resolution optical imaging system, 
by which single atoms are detected with near-unity fi- 
delity on individual sites of an OL [53[. A similar mea- 
surement has been performed for revealing phase transi- 
tions of an atomic Ising chain [54j. 

In summary, we have studied the superfluidity and 
magnetic properties of ultracold bosons with the syn- 
thetic 3D SO coupling in a square OL. The lowest en- 
ergy band exhibits one, two or four pairs of degener- 
ate single-particle ground states depending on the SO- 
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coupling strengths, which can lead to the condensate 
states with spin-stripes for the weak atomic interactions. 
In the MI state with one particle per site, an effective spin 
Hamiltonian with the 3D DM spin interactions is derived. 
The spin Hamiltonian with an additional Zeeman field 
has a rich phase diagram with spiral, stripe, vortex crys- 
tal, and Skyrmion crystal spin-textures in each xy-plane 
layer. The Skyrmion crystals extended along the z axis 
in a columnar manner can be tunable with square and 
hexagonal symmetries, and stable against the inter-layer 
spin-spin interactions in a large parameter region. 
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